#主要实现：
#获得主城区范围
#不同类型城市，聚类参数可调
#生成包含主城区的最小凸多边形

import os
import utm
import pyproj
import pandas as pd
from osgeo import gdal
from geopandas import GeoDataFrame
from shapely.geometry import Point
from shapely.geometry import MultiPoint
from shapely.ops import transform

#栅格处理函数，输入矩形经纬度坐标
# 返回zone_point和pop_border(区域的几何点和人口降序排列的15%点的值)
def raster_operation(xmin, ymax, xmax, ymin,filename):
    # 读取部分数据集
    print('******************')
    print('开始读取区域栅格...')
    # os.chdir(r'E:\MyWork\0_SD_Local\project\03_cq\data')
    in_ds = gdal.Open(filename)
    in_gt = in_ds.GetGeoTransform()  # 获得输入栅格数据的地理变换（geotransform）

    inv_gt = gdal.InvGeoTransform(in_gt)  # 获得地理变换（geotransform)的逆
    if gdal.VersionInfo()[0] == '1':
        if inv_gt[0] == 1:
            inv_gt = inv_gt[1]
        else:
            raise RuntimeError('Inverse georansform failed')
    elif inv_gt is None:
        raise RuntimeError('Inverse geotransform failed')

    offsets_ul = gdal.ApplyGeoTransform(inv_gt, xmin, ymax)  # 获得左上角和右下角偏移坐标（图像坐标）
    offsets_lr = gdal.ApplyGeoTransform(inv_gt, xmax, ymin)
    off_ulx, off_uly = map(int, offsets_ul)
    off_lrx, off_lry = map(int, offsets_lr)
    print('左上角图像坐标:')
    print(off_ulx, off_uly)
    print('右下角图像坐标:')
    print(off_lrx, off_lry)

    rows = off_lry - (off_uly - 1)  # 获得要提取区域的行数和列数
    colums = off_lrx - (off_ulx - 1)  # 图上坐标中左上角的x,y坐标最小
    print('行数：' + str(rows), '列数：' + str(colums))

    in_band = in_ds.GetRasterBand(1)
    data = in_band.ReadAsArray(off_ulx, off_uly, colums, rows)  # 研究区域的栅格数据，Numpy二维数组格式
    # print("区域栅格数组：\n"+str(data))

    # 研究区域的栅格数据处理
    print('开始栅格数据处理...')
    xyCoors = []
    GRID_CODE = []
    for r in range(off_uly, off_lry + 1):  # 获得区域栅格的中心点坐标以对应的值
        for c in range(off_ulx, off_lrx + 1):
            # xyCoors.append(src.xy(r, c))
            xyCoors.append(gdal.ApplyGeoTransform(in_gt, c, r))
            GRID_CODE.append(data[r - off_uly, c - off_ulx])

    print('坐标点个数：')
    print(len(xyCoors))
    print('点的值个数：')
    print(len(GRID_CODE))

    s1 = pd.Series(xyCoors)  # 创建DataFram
    s1.name = 'geometry'
    s2 = pd.Series(GRID_CODE)
    s2.name = 'GRID_CODE'
    df = pd.DataFrame([s1, s2], index=['geometry', 'GRID_CODE'])
    df = df.T

    geometry = [Point(xy) for xy in df['geometry']]  # DataFram转化为GeoDataFram
    df = df.drop(['geometry'], axis=1)  # 从行或列中删除指定的标签。
    crs = {'init': 'epsg:4326'}
    zone_point = GeoDataFrame(df, crs=crs, geometry=geometry)

    L = sorted(GRID_CODE,reverse=True)
    n = int(len(xyCoors)*0.15)
    pop_border = L[n]
    print('栅格处理已完成!')
    print('******************')
    return zone_point,pop_border

#文件夹目录处理
china_population_2016_Dic = "../files/china_population_2016.tif"

f = open("../files/runinfo_vc6.txt")
next(f)
for d in f:
    d = d
coorDic = str(d)+"/data/coor.txt"
latlon_utmDic = str(d)+"/data/latlon-utm.txt"
mainAreaDic = str(d)+"/cache/mainArea.txt"
autoZoomParaDic = str(d)+"/cache/autoZoomPara.txt"

#灯光人口数据文件路径
f = open(autoZoomParaDic)
next(f)
d2 = f.readline()
d3 = f.readline()
d4 = f.readline()
pop_lig_type = int(d2.split()[1])
pop_lig_file = d3.split()[1]
city_type = d4.split()[1]

#获得区域边界矩形的左上角和右下角经纬度坐标
coor_y,coor_x = [],[]
lat,lon = [],[]
for line in open(coorDic):                        #coor.txt中点的编号不连续，总共有6425个
    coor_y.append(float(line.split()[1]))
    coor_x.append(float(line.split()[2]))

if(len(coor_y) == len(coor_x)):                                                                     # 获得点的个数
    codeNum = len(coor_y)
else:
    print("coor error")

f = open(latlon_utmDic)                          #获得zone和band
d1 = f.readline()
zone_number = int(d1.split()[5])
zone_letter = d1.split()[6]

for i in range(0,codeNum):                              #UTM转经纬度
    lat_temp, lon_temp = utm.to_latlon(coor_y[i], coor_x[i], zone_number, zone_letter)
    lat.append(lat_temp)
    lon.append(lon_temp)
xmin,ymax = min(lon),max(lat)                           #矩形左上角经纬度坐标
xmax,ymin = max(lon),min(lat)                           #矩形右下角经纬度坐标
print('研究区域矩形范围：')
print('左上角坐标:')
print(xmin,ymax)
print('右下角坐标:')
print(xmax,ymin)

#调用栅格数据处理函数
#判断输入的是灯光还是人口
if (pop_lig_type == 0):
    zone_point,pop_border = raster_operation(xmin, ymax,xmax, ymin,pop_lig_file)
    print("输入的是人口数据")
else:
    zone_point,pop_border = raster_operation(xmin, ymax,xmax, ymin,china_population_2016_Dic)
    print("输入的是灯光数据")

print("zone_point已生成!")
# print(zone_point)

#人口平均值作为阈值筛选点
print('开始阈值筛选...')
print('人口阈值：'+str(pop_border))
zone_point_border = zone_point[zone_point.GRID_CODE.values>pop_border]
print('zone_point_border已生成！')
# print(zone_point_border)

dataSet = []
lon_border = []
lat_border = []
for point in zone_point_border.geometry:
    list(point.coords)
    lon_border.append(point.x)
    lat_border.append(point.y)
    dataSet.append([point.x,point.y])
print('zone_point_border坐标信息已提取！')

#聚类分析
print('开始聚类分析...')
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import DBSCAN
from math import *
from scipy.spatial.distance import squareform,pdist
from collections import Counter

#加载数据集
X = np.asarray(list(dataSet))
print('坐标点的个数：'+str(len(X)))

#根据两点的经纬度计算实际距离，单位：千米
def haversine(lonlat1, lonlat2):
    lat1, lon1 = lonlat1
    lat2, lon2 = lonlat2
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])

    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2
    c = 2 * asin(sqrt(a))
    r = 6371                                                                  #地球半径
    return c * r

distance_matrix = squareform(pdist(X, (lambda u, v: haversine(u, v))))       #两点之间距离矩阵
d = np.array(distance_matrix)
print('距离矩阵型号:'+str(d.shape))

#聚类分析
i = int(city_type)                    #获取用户选择的城市类型,int  特大：0、大：1、中：2、小：3
epsList = [1.414,1.414,1,1]
samplesList = [8,6,4,2]
db = DBSCAN(eps = epsList[i],min_samples = samplesList[i],metric='precomputed').fit(d)      #默认值eps = 0.5,min_samples = 5. #10.8  7
core_samples_mask = np.zeros_like(db.labels_, dtype=bool)              #某点的2公里范围内有5个点，则此点作为核心点
core_samples_mask[db.core_sample_indices_] = True
labels = db.labels_
n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)
print('聚类已完成！')
print('标签个数: %d' % len(labels))
print('类的个数: %d' % n_clusters_)

# #聚类前原图
# plt.scatter(X[:, 0], X[:, 1], marker='o')
# plt.show()
#
# # 聚类结果画图
# unique_labels = set(labels)
# colors = [plt.cm.Spectral(each)
#           for each in np.linspace(0, 1, len(unique_labels))]
# for k, col in zip(unique_labels, colors):
#     if k == -1:
#         # Black used for noise.
#         col = [0, 0, 0, 1]
#     class_member_mask = (labels == k)
#     xy = X[class_member_mask & core_samples_mask]
#     plt.plot(xy[:, 0], xy[:, 1], 'o', markerfacecolor=tuple(col),
#              markeredgecolor='k', markersize=8)
#
#     xy = X[class_member_mask & ~core_samples_mask]
#     plt.plot(xy[:, 0], xy[:, 1], 'o', markerfacecolor=tuple(col),
#              markeredgecolor='k', markersize=6)
#
# plt.title('Estimated number of clusters: %d' % n_clusters_)
# plt.show()

#找出核心点最多的类别中的点
result = Counter(labels)                                                 #统计各类中点的数量
result[-1] = 0                                                           #异常点的数量设置为0，防止其数量本来最大
print(result)
mostLabel = max(result,key=result.get)                                   #找出核心点最多的类中点的labels


name = list(labels)                                                      #查找list中的某个元素的所有的下标
first_pos = 0
poiLocation = []
for i in range(name.count(mostLabel)):                                  #统计largestCluster出现的次数作为范围
    new_list = name[first_pos:]
    poiLocaTemp = first_pos + new_list.index(mostLabel)
    poiLocation.append(poiLocaTemp)                                     #poiLocation存储mostLabel在列表中的位置
    next_pos = new_list.index(mostLabel) + 1
    first_pos += next_pos
print('mostLabel下标个数：'+str(len(poiLocation)))

new_dataSet = []                                                        #new_dataSet存储one_point_max经纬度坐标
for p in poiLocation:
    new_dataSet.append(dataSet[p])
print('zone_point_max的个数:'+str(len(new_dataSet)))

#经纬度转UTM
def wgs2proj(y_84,x_84):
    (x_coor, y_coor, temp_1, temp_2) = utm.from_latlon(x_84, y_84)                  #此处的x_coor和y_coor应该调换
    if not (std_temp_1 == temp_1):                                                 # 存在跨zone的点
        p1 = pyproj.Proj(proj='utm', zone=temp_1, datum='WGS84')
        p2 = pyproj.Proj(proj="utm", zone=std_temp_1, datum='WGS84')
        x_coor, y_coor = pyproj.transform(p1, p2, x_coor, y_coor)                  #转换带号
    return x_coor, y_coor

#确定包含主城区的最小凸多边形
std_temp_1 = zone_number
print("--开始确定主城区范围--")

# print(new_dataSet)
m = MultiPoint(new_dataSet)
mainAreaPoly = m.convex_hull
# print(mainAreaPoly)

#transform（func，geom ）
# 将func应用于geom的所有坐标，并从变换的坐标返回相同类型的新几何
mainAreaPoly_utm = transform(lambda x, y: wgs2proj(x, y), mainAreaPoly)
# print(mainAreaPoly_utm)
mainArea = list(mainAreaPoly_utm.exterior.coords)#凸多边形边界坐标
# print(mainArea)
XY = np.asarray(mainArea)
XY = XY*[1000]

with open(mainAreaDic, 'w') as f:
    f.write('边界点x坐标 边界点y坐标\n')
    for i in range(len(XY)):
        for j in range(len(XY[i])):
            f.write(str(XY[i][j]))
            f.write(' ')
        f.write('\n')
f.close()
print("--包含主城区的凸多边形已生成--")


